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ABSTRACT 


In  this  study,  several  turbulence  models  are  analyzed 
in  a  3-D  finite  element  model  of  a  micro-heat  exchanger. 
The  micro-heat  exchanger  consists  of  a  narrow  planar  flow 
passage  between  parallel  walls  with  small  cylindrical  pin 
fins  spanning  these  walls  with  axes  perpendicular  to  the 
direction  of  flow.  Turbulence  model  performance  is 
compared  with  baseline  experimental  data  available  in  the 
literature  that  cover  a  range  of  low  turbulent  Reynolds 
numbers  and  spacing  configurations.  The  metric  for  these 
comparisons  is  an  array  averaged  Nusselt  Number. 
Adjustments  made  to  the  coefficients  in  the  turbulence 
models  are  explained  in  terms  of  their  physical 
significance  to  the  complex  flow  environment  of  a  pin  fin, 
cross  flow,  micro-heat  exchanger.  Applications  of  this 
research  include  cooling  of  turbine  blades  and  of  closely 
spaced  electronics. 
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I. 


INTRODUCTION 


High  temperatures  on  power  turbine  blades  have  been  a 
limiting  factor  in  advancing  the  design  of  gas  turbine 
engines.  Conventional  cooling  schemes  in  advanced  gas 
turbine  engines  will  be  unable  to  remove  the  heat 

transferred  from  the  mass  stream  fast  enough  to  prevent 

damage  to  the  blades.  A  new  cooling  scheme  has  been 
proposed  which  may  be  adaptable  to  the  turbine  blade 

cooling  problem.  It  consists  of  a  narrow  flow  passage 

between  two  endwalls  is  separated  by  rows  of  short  pin  fins 
with  axes  perpendicular  to  the  coolant  flow.  It  is 
envisioned  that  this  flow  passage  would  be  tightly 

"wrapped"  around  a  turbine  blade  forming  a  protective 
cooling  shroud.  Bleed-off  air  for  cooling  would  be 
supplied  through  the  blade  root  and  distributed  at  the 

leading  edge  of  the  passage  flow  with  an  exit  to  the 

turbine  flow  at  the  trailing  edge  of  the  blade.  This  paper 
addresses  the  turbulence  models  used  in  the  3-D  numerical 
modeling  of  two  arrangements  of  this  cooling  scheme  over 
several  low,  but  fully  turbulent  Reynolds  numbers.  The 

turbulence  models  have  been  fine  tuned  to  agree  with 
available  experimental  data  in  the  literature,  and  optimal 
models  and  coefficients  have  been  identified. 
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II.  BACKGROUND  AND  OBJECTIVES 


A.  HEAT  EXCHANGER  CONSTRUCTION 

The  micro-heat  exchanger  modeled  in  this  study  is  made 
up  of  a  narrow  flow  passage  supported  by  10  staggered  rows 
of  pin  fins  perpendicular  to  the  flow.  Figure  1  is  a 
definition  sketch  of  a  typical  configuration,  similar  to 
the  type  constructed  by  Metzger,  et  al.  (1982)  which  forms 
the  basis  of  corroboration  for  the  current  numerical  study. 
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Figure  1 .  Definition  sketch  of  a  staggered  pin- 

fin  array 

The  pins  are  arranged  with  a  constant  axial  pitch  (X, 
distance  between  pin  centers  in  the  flow  direction)  to  form 
the  reference  length  scale.  The  axial  pitch  was  chosen  to 
be  12.7mm  to  be  consistent  with  earlier  studies.  Two 
different  diameter  pins  were  used  in  this  setup,  8.46mm  and 
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5.08mm.  The  X/D  and  S/D  ratios  define  the  relative 
streamwise  and  spanwise  pin  spacing  respectively  for  a 
general  model.  As  a  result  of  these  settings  the 

transverse  pitch  (S)  was  not  a  constant  from  model  to 
model,  although  axial  pitch  was. 

As  in  the  work  of  Metzger  et  al  (1982),  the  endwalls 
(the  broad  upper  and  lower  plates)  were  modeled  as 
isothermal  surfaces.  The  pins  were  modeled  as  isothermal 
at  the  temperature  of  the  bounding  endwalls  as  a  first 
approximation.  It  can  be  shown  with  a  standard  pin  fin 
relation  (Incropera  and  Dewitt,  1996)  that  the  isothermal 
assumption  for  the  pins  is  a  valid  assumption.  The  side 
walls  were  modeled  as  insulated,  again  to  duplicate 
Metzger's  work  and  an  insulated  inlet  and  outlet  section 
were  added  (not  shown  in  Figure  1) . 

B .  PREVIOUS  WORK 

The  staggered  pin-fin  array  geometry,  in  one  form  or 
other,  serves  as  the  basis  for  a  large  number  of  compact 
heat  exchanger  configurations,  and  has  received 
considerable  attention  in  the  past.  The  experimental  work 
of  VanFossen  (1982),  Metzger  et  al  (1982),  Chyu  (1990), 
Chyu  and  Goldstein  (1991),  and  Jubran  et  al .  (1993)  and 

others  have  explored  various  issues  such  as  the  differences 
between  short  and  long  pins,  the  contribution  of  pin  heat 
transfer  versus  endwall  heat  transfer,  the  presence  of  an 
optimal  streamwise  and  spanwise  spacing/pitch,  etc. 

Only  recently,  with  the  advent  of  more  powerful 
computational  machines  and  software  tools,  have  numerical 
studies  of  compact  heat  exchangers  been  attempted  as 
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reviewed  by  Shah  et  al.  (2001) .  For  the  staggered  pin  fin 
array  configuration  in  particular,  Donahoo  et  al .  (2001) 
have  used  a  general  purpose  viscous  flow  solver  to  simulate 
the  fluid  flow  and  heat  transfer  behavior  in  a  2-D 
numerical  study.  However  the  2-D  nature  cannot  be  used  to 
treat  the  endwalls  where  the  pins  meet  the  surface. 

The  experimental  baseline  for  this  study  is  the 
Metzger  et  al .  work  of  1982.  Metzger  showed  the  variation 
of  streamwise  heat  transfer,  overall  array  heat  transfer, 
and  overall  flow  friction  in  staggered,  short  pin  fin 
arrays.  His  experimental  work  has  provided  the  baseline 
for  this  computer  model.  His  experimental  work  showed  that 
there  was  a  substantial  difference  between  short  pin  fin 
arrays  and  long  cylinder  (tube  bank)  arrays.  Long 
cylinders  had  been  the  previous  model  for  the  array.  The 
relations  Metzger  developed  for  overall  heat  transfer 
showed  a  lower  array  averaged  Nusselt  number  over  the  array 
based  on  his  experiments.  The  long  cylinder  predictions 
did  not  correctly  predict  the  convection  coefficients 
observed  in  the  short  pin  fin  array. 

C .  OBJECTIVE 

The  primary  objective  of  this  study  was  to  select  a 
turbulence  model  that  accurately  reflects  the  experimental 
results  seen  in  previous  work.  As  part  of  the  selection 
process,  the  turbulence  model  had  to  be  optimized  since  the 
default  coefficients  in  a  model  are  typically  better  suited 
for  flat  plate  or  long  cylinder  tube-bank  type  open 
geometries,  but  not  very  applicable  to  confined  tortuous 
flow  geometries  such  as  the  short  pin  fin  array  in  this 
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study.  Another  objective  is  to  leave  the  code  open  for 
future  work  where  pin  shape  and  spacing/layout  could  be 
changed.  The  goal  would  be  to  arrive  at  a  recommendation 
for  the  most  optimal  turbulence  model/s  and  their 
coefficients  that  would  provide  the  best  corroboration  with 
available  experimental  data.  Furthermore  it  is  hoped  that 
subject  to  future  experimental  work,  the  current  findings 
can  be  extended  to  have  validity  over  a  wide  range  of 
Reynolds  numbers  and  geometrical  configurations  thus 
providing  greater  applicability  and  confidence  in  future 
numerical  simulations. 
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TURBULENCE  MODELS 


A.  INITIAL  MODELING 

1 .  Boundary  Conditions 

The  numerical  model  of  Metzger's  experiment  consists 
of  two  different  boundaries.  Fluid  elements  were  defined 
which  had  rigid  boundaries,  eliminating  the  need  for 

defining  a  rigid  structure  to  surround  the  fluid  elements. 
To  reduce  the  number  of  nodes  required  for  modeling, 

symmetry  was  exploited.  A  horizontal  plane,  parallel  to 
the  flow  direction,  halfway  between  the  upper  and  lower 
isothermal  boundaries  of  the  passage  defines  an  adiabatic 
boundary.  No  net  properties  cross  this  boundary.  The 
length  scale  in  the  flow  passage  is  sufficiently  small  that 
buoyancy  effects  are  negligible.  This  boundary  was  defined 
as  insulated  and  velocity  in  the  vertical  direction  was 

zero.  However  the  no-slip  condition  along  the  plane  was 
released.  A  similar  symmetry  plane  was  defined 

perpendicular  to  the  first  symmetry  plane  and  parallel  to 
the  flow  direction.  The  boundary  was  again  modeled  as 
insulated,  crossing  velocity  was  set  to  zero  and  the  no¬ 
slip  condition  along  the  plane  was  released.  These  two 
planes  allowed  a  model  to  only  have  nodes  in  one  quadrant 
as  the  other  quadrants  would  produce  the  same  result. 
Observing  the  downstream  direction,  the  lower-right 
quadrant  was  the  quadrant  chosen  for  modeling. 
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Symmetry 


Figure  2.  Modeled  Portion  of  a  Pin  Fin  Heat 

Exchanger 
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Figure  3.  Actual  model  used  for  experiment 

The  pins  were  modeled  as  isothermal  and  so  was  the 
supporting  endwall.  The  remaining  wall,  the  sidewall, 
parallel  to  the  pin  axes  was  insulated.  The  no  slip 
condition  was  applied  to  all  boundaries  except  the  symmetry 
planes.  Straight  duct  entrance  and  exit  regions  were 
attached  at  the  inlet  and  outlet  sections  of  the  array  test 
section  to  be  consistent  with  the  experimental  rig  of 
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Metzger  et  al  (1982)  .  They  were  modeled  slightly  longer 
than  the  Metzger  model  to  prevent  the  outlet  conditions 
from  affecting  the  test  section. 

The  isothermal  boundaries  were  fixed  at  306  K.  Dry 
air  (Pr  =  0.7)  was  used  as  working  fluid,  entering  the 

inlet  section  at  300K,  p  =1.1769-^-,  |l  =  1.8468xl0”5 — — , 

m  m  ■  s 

.  w  j 

K  —  2.6294x10  " - ,  and  C„  =1004.0 - .  Inlet  velocity  was 

m-K  p  kg-K 

specified  at  the  entrance  to  the  test  section  based  on  a 
maximum  velocity  type  Reynolds  number. 

The  maximum  velocity  Reynolds  number  is  defined  as 

Re  =  ^ D  (Metzger  et  al .  ,  1982),  where  Vmax  is  the 

(i 

inviscid,  geometry  based  maximum  velocity  the  fluid  reaches 
when  passing  through  the  test  section.  D  is  the  pin 
diameter . 

2 .  Initial  Conditions 

Each  run  was  initiated  with  the  walls  and  pins  of  the 
test  section  at  306K  and  the  incoming  air  at  300K.  Inlet 
velocity  was  determined  by  the  specified  Reynolds  number 
and  fluid  properties. 

3 .  Mesh  Structure 

The  mesh  was  specified  to  concentrate  nodes  in  the  no 
slip  boundary  areas  of  the  model,  where  flow  and 
temperature  gradients  might  be  expected  to  be  highest  (Fig 
4)  .  Around  the  pins,  nodes  were  concentrated  in  the  high 
velocity  regions  where  conditions  change  the  most  (Fig  5)  . 
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Se='?310,S/D=2.5,m=l.S,TW  =306 ,  h/v/p=5/16/4(QXG.250,SKB 


Figure  4 .  Standard  Mesh  showing  nodes 

concentrated  in  the  areas  of  the  no  slip 

boundary . 


Figure  5.  Standard  Mesh  (Top  View)  with  nodes 

concentrated  around  the  high  velocity  regions  of 

the  pins 
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Nodes  were  also  concentrated  along  the  endwall. 


For  constructing  the  numerical  model  of  the  heat 
exchanger  and  for  the  majority  of  the  data  gathering,  the 
finite  element  based  modeling  and  simulation  package,  ANSYS 
was  used.  The  versions  available  during  this  research  were 
ANSYS  version  5.7,  6.0,  and  6.1.  The  bulk  of  the  research 
was  done  using  ANSYS  6.0.  The  ANSYS  license  at  the  Naval 
Postgraduate  School  (NPS)  limited  models  to  256,000  nodes. 
Several  runs  were  attempted  using  the  Department  of  Defense 
High  Performance  Computing  Major  Shared  Resource  Center 
(HPCMSRC) .  The  ANSYS  license  at  the  HPCMSRC  was  node 
unlimited. 

At  the  entrance  and  exit  of  the  test  section,  as  well 
as  midway  between  each  row  of  pins  (apparent  in  figures  4 
and  5),  a  vertical  plane  of  nodes  was  constructed.  This 
allowed  bulk  temperature  and  mass  flows  to  be  calculated 
independent  of  the  ANSYS  program.  Such  a  layer  of  nodes 
would  also  facilitate  future  work  where  pin  shape  and 
layout  are  expected  to  be  changed. 


B .  SKE  MODEL 

The  Standard  k-8  (SKE)  Model  is  the  default  turbulence 
model  for  ANSYS  6.0.  The  SKE  model  is  a  Reynolds  Averaged 
Navier-Stokes (RANS)  type  of  Eddy  Viscosity  Model (EVM) .  The 
SKE  model  uses  2  equations  to  close  the  governing 
equations.  The  EVM  approach  introduces  a  turbulent 
viscosity  term,  (I, ,  which  relates  the  fluctuating  velocities 
of  a  turbulent  flow  to  a  viscous  stress  term. 
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in  tensor  notation,  which  is  used  in  the 


—UjU 


i 


(I,  dui 

p  dXj  ' 


Reynolds  stress  terms  of  the  time  averaged,  turbulent, 
incompressible  Navier-Stokes  equations,  as  well  as  in  the 
energy  equation  and  the  continuity  equation.  As  an 
example,  the  x  component  of  the  Reynolds  stress  term  is: 


<b  )-^(pft)~p(pVX) '  "here  APf  are  the 

fluctuating  components  of  the  velocities.  The  Reynolds 
stress  terms  contain  the  indeterminate  part  of  the  Navier- 
Stokes  relations,  the  mean  of  the  product  of  the 
fluctuating  velocity  components.  Replacing  that  product 
with  the  relation  involving  the  turbulent  viscosity  allows 
the  Navier-Stokes  momentum  equations,  continuity  and  energy 
equations  to  be  reshaped.  In  the  SKE  model  two  transport 
equations  are  derived  from  the  governing  equations,  the 
Turbulent  Kinetic  Energy  equation: 

apt  | a(pv^)  fdPtd+hPth 

dt  dx  dy  dz 


d_ 

dx 


-pe  + 


|I,  dk 

c)x  j 

QP  n, 


+  ■ 


dy 


(I,  dk 
Gk  dy 


+ 


dz 


dT 


dT 


(i,  dk 

SzJ 
dT  A 


+  M> 


(1.1) 


and  the  Dissipation  Rate  equation: 


3pe  |  3(pVe)  |  ^(pKe)  |  3(pKe) 

dt  dx  dy  dz 


d 

'  (i,  3e  ^ 

d 

J _ 

V,  de  ^ 

d 

_| _ 

f  (i,  de  N 

dx 

aP  dx 

K  £  J 

dy 

Of  dy 

K  £  ■  ) 

dz 

G„  dz 

v  £  y 

+  C^tj<S>-C2p^- 
k  k 
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where  k  is  the  turbulent  kinetic  energy,  Vx,Vy,Vz  are 
the  instantaneous  (non-fluctuating)  velocities,  o k  is  the 

Schmidt  number  for  the  transport  of  turbulent  kinetic 
energy,  $  is  the  viscous  dissipation  term 

f du ,  dut  ^  .  ,  ,  .  . 

<J>  =  pi  — — -  ,  e  is  the  turbulent  kinetic  energy 

[dxk  Pxk  y 

dissipation  rate,  C3  and  C4  are  the  buoyancy  constants,  P 
is  the  coefficient  of  thermal  expansion,  a,  is  the  turbulent 
Prandtl  number,  g  is  the  acceleration  scalar,  T  is  the 
static  temperature,  Oe  is  the  Schmidt  number  for  the 
transport  of  turbulent  kinetic  energy  dissipation  rate,  ClE 
is  the  shear  rate  generation  term  multiplier,  C 0  is  scaled 
diffusivity  rate  multiplier,  C is  the  turbulent  viscosity 
constant,  and  (I,  is  the  turbulent  viscosity  broken  down  in 

p  k~ 

the  equation  to  (I,  =  C  .  Table  1  contains  the  values  of 

£ 

the  SKE  constants: 
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1  .  0 

Q 

1  .  0 

C4 

0  .  0 

P 

0  .  0 

1  .  0 

°£ 

1 . 3 

Qe 

1.44 

1 . 92 

C, 

0.09 

Table  1.  Constants  for  SKE  Model 


Solving  these  equations  at  each  node  in  the  model  for 
a  set  number  of  iterations  yields  values  for  k  and  £  which 
are  used  to  calculate  |1,  from  the  previous  equation.  (lf  is 

used  to  calculate  the  indeterminate  fluctuating  velocities 
in  the  Navier-Stokes  equations.  The  Navier-Stokes 
equations  are  then  solved  for  a  new  set  of  mean  velocities. 

This  new  set  of  velocities,  with  fluctuating 
components,  is  then  used  to  solve  for  temperature  in  a 
first  law  relation: 


pc 


( dT  dT  dT  dT  ) 

—  +  y  — +  V  — +  V  — 

dt  x  dx  y  dy  z  dz 


=  K 


d2T 
dx 7 


-  + 


v 


d2T  d2T ' 

dy2  dz2 


(1.3) 
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where  K  is  the  thermal  conductivity  of  the  fluid, 
assumed  to  be  isotropic,  c  is  the  specific  heat,  and  the 
velocities  V  ,V  ,V  contain  the  mean  and  fluctuating 

X7  y7  Z 

components . 

The  local  heat  flux,  q  ,  can  be  estimated  at  non- 

insulated  or  adiabatic  boundaries  by  the  following  relation 

where  y  is  the  direction  normal  to  the  transfer  surface. 

From  this,  Newton's  Law  of  Cooling  leads  to  a 

y=0 

convection  coefficient,  h  =  ~, — — - -r ,  where  T„  is  the  surface 

(f-t) 

temperature  and  is  the  free  stream  temperature.  Since 

the  flow  complex,  ANSYS  6.0  uses  TB  in  place  of  7^.  TB  is 

the  bulk  temperature  of  the  adjacent  fluid  element. 

The  ultimate  goal  of  the  simulation  is  to  calculate  an 
array  averaged  Nusselt  number,  so  an  overall  array 
convection  coefficient,  h ,  is  calculated  based  on  the  total 
heat  transferred  to  the  fluid. 


vdT\ 
q  =  —  K  — 

ay1 


h  = - - -  (1.4) 

A  ■  LMTD 

where  q  is  the  total  amount  of  heat  in  watts 
transferred  across  the  boundaries,  A  is  the  wetted  area  of 
the  test  section,  and  LMTD  is  the  Log-Mean  Temperature 
Difference . 


LMTD  = 


In 


f 


AT 

T  -T. 

wall _ m_ 

T  -T 

wall  out 


( 1.5 ) 
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is  the  bulk  inlet  temperature 


AT  is 


where 


T 


calculated  from  AT  — - ,  T  ,,  is  the  constant  isothermal 

•  /-i  wall 

m  •  C 

p 

wall  temperature  for  both  endwalls  and  pins,  and 
T  -T  +  AT  .  The  array  averaged  Nu  is  calculated  from: 

out  in  -1  -1 


(1.6) 


where  D  is  the  diameter  of  a  pin  and  K  is  the  thermal 

T  +T 

conductivity  of  air  evaluated  at  — - n^!L  m  This  formulation 

2 

of  Nusselt  number  is  the  most  appropriate  for  comparing  the 
experimental  results  of  Metzger  with  the  current  ANSYS 
simulations . 

One  drawback  to  the  SKE  model,  as  well  as  all  the  2 
equation  EVM's,  is  that  it  does  not  accurately  predict 
turbulent  kinetic  energy  and  dissipation  rate  at  elements 
along  the  wall.  In  this  case  a  wall  turbulence  model  is 
used.  The  ANSYS  default  model  for  wall  turbulence  is  the 
Van  Driest  model.  The  Van  Driest  model  is  used  at  the  wall 
nodes  to  determine  viscosity.  It  assumes  an  eddy  viscosity 
of  zero  at  the  wall  and  gives  an  accurate  value  for  wall 
conduct ivity .  ANSYS  6.0  uses  the  continuous  Van  Driest 
formulation  with  the  velocity  at  a  certain  distance  from 
the  wall  as  a  reference  velocity.  This  form  of  the  log  law 
of  the  wall  is  solved  for  a  wall  viscosity  from 


where  X  is  the  shear  stress,  K  is  the  slope  of  the 
wall  parameter,  E  is  the  Law  of  the  Wall  constant,  |i,  is 

the  laminar  absolute  viscosity,  and  5  is  the  distance  from 
the  wall.  The  wall  element  effective  viscosity,  |iu, ,  is 
derived  from: 


tan 


(1.8) . 


C.  RENORMALIZED  GROUP  (RNG)  MODEL  -  DIFFERENCES  FROM  SKE 

The  RNG  model  uses  an  advanced  statistical  technique 
called  Renormalization  Group  Theory  applied  to  the 
instantaneous  Navier-Stokes  equations.  The  two  transport 
equations  derived  for  the  SKE  model  appear  in  the  RNG  model 
with  the  exception  of  the  ClE  multiplier  for  the  kinetic 
energy  generation  term  in  the  dissipation  rate  equation. 
In  the  RNG  model  the  ClE  term  takes  the  form 


C,  =1.42- 

le 


*1 


_JL 


l+|3rr 


(1.9) 


where  ^=4.38,  (3=0.012 

stress  deformation  tensor, 
function  of  velocities  in  the 


S  =  — 

iJ  2 


dV  dV 
— —+ — - 
dV  dV 

J  '  ) 


(1.10) 


,  and  T|  is 
S...  The 

ij 

f lowf ield. 


a  function 
tensor,  S  , 


of  the 
is  a 


There  are  differences  in  the  buoyancy  terms  of  the 
equations  but  for  this  heat  exchanger  as  noted  earlier 
buoyancy  effects  can  be  ignored.  In  the  RNG  model  several 
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constants  have  different  numerical  values 


the  default 


numerical  constants  are  summarized  in  Table  2: 


0.72 

Q 

1  .  0 

Q 

o 

o 

P 

0 .012 

1 . 0 

0.72 

TL 

4.38 

1 . 68 

c* 

0 .085 

Table  2.  RNG  Model  Coefficients 

This  gives  the  RNG  model  an  added  level  of  complexity 
and  capability  over  the  SKE  model  due  to  the  more 
distributed  nature  of  accounting  for  the  turbulent  kinetic 
energy  generation  term  in  the  flow  field.  In  the  SKE 
model,  the  kinetic  energy  generation  coefficient  in  the 
dissipation  rate  equation  is  the  same  throughout  the 
flowfield.  In  the  RNG  model,  the  coefficient  is  dependent 
upon  the  velocity  components  at  each  point  in  the  flow 
field.  ANSYS  6.0  uses  a  slightly  more  simplified  version 
than  that  presented  in  Numerical  Analysis  of  Compact  Heat 
Exchanger  Surfaces  (R.  K.  Shah  et  al .  ,  2001)  . 
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Past  research,  as  noted  R.  K.  Shah  et  al .  ,  describes 
RNG  as  being  more  accurate  and  reliable  than  other  2 
equation  EVM  models.  Because  it  is  a  new  model,  relatively 
few  studies  have  been  done  with  it,  but  it  appears  to  do 
better  than  SKE  with  low  Reynolds  numbers  but  worse  when 
vortex  shedding  becomes  an  issue,  such  as  with  bluff  bodies 
or  flow  around  cylinders  as  in  external  flows  with  open 
geometries.  Although  the  comparisons  were  based  on  2D 
models  only,  RNG's  success  appears  to  be  highly  geometry 
dependent  . 

D.  NKE  -  DIFFERENCES  FROM  SKE 

The  New  k  —  £  Model  (NKE) ,  also  known  as  the  Realizable 
k  —  8  Model,  was  proposed  by  T.  H.  Shih  et  al .  in  1994.  The 
model  was  primarily  designed  for  high  Reynolds  number  flows 
but  has  had  some  success  for  lower  Reynolds  number  flows  as 
well.  The  significant  differences  for  a  micro  heat 
exchanger  type  of  flow  environment  are  that  the  NKE  model 
uses  a  different  formulation  for  the  C,,  constant  and  a 

different  form  of  the  kinetic  energy  dissipation  rate 
equation . 

The  turbulent  viscosity  (I,  is  still  calculated  the  same 
way  but  the  C  term  is  now  a  function  of  deformation 
tensors  and  antisymmetric  rotation  tensors.  ANSYS  6.0 
replaces  C with  this  function: 


where  the  C  Cl  8  term  involves  the  angular  velocity  of  the 

r  m  mij 

coordinate  system  which  for  this  application  is  not 
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rotating,  so  that  term  drops  out.  The  entire  term  under 
the  radical  is  the  mean  velocity  in  the  law  of  the  wall 
model . 


The  dissipation  rate  equation  now  becomes: 

3pe  |  a(pV£)  |  a(pyve)  |  d(pKe) 

dt  dx  dy  o k 


d_ 

dx 


}L^d£_ 

v°e  dx 


d 


dy 


(i,  de 
°E  dy 


+ 


dz 


f  \ i,  3e  ^ 
30  dz 


+  pclte  j_. 


3V,  „  ' 
3  V  + 

V  1  J 


(dyi  ..  A 

dsr 

]  j 


£ 2  c  (i-c3)ppfcf  ar  ar  ar^ 

ar  ay  az 


-C2p^-+  * 
k 


(1.12) 


and  the  C,  multiplier  now  becomes 


C,  =  max 

le 


f  r  ^  ^ 

v  1M  'tl  +  5  j 


(1.13) 


where , 


and  C 


dV  dV . 


•+- 


ay.  ay  , 

'  '  j 


dv,  avp 

ay.  ay.  , 

.  1  1 ) 


is  a  constant 


(1.14) 


The  numerical  values  for  the  default  NKE  model  are 
summarized  in  Table  3: 


°A 

1  .  0 

Q 

1  .  0 

Q 

o 

o 

P 

o 

o 
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1  .  0 

°£ 

1.2 

1 . 90 

C\M 

0.43 

Table  3.  NKE  Model  Coefficients 

E.  FAILURES  -  ZEM,  GIR,  SZL 

There  were  three  other  turbulence  models  that  were 
evaluated  in  this  study,  a  zero  equation  model  (ZEM)  ,  the 
non-linear  model  of  Girimaji  (GIR) ,  and  the  Shih,  Zhu, 
Lumley  model  (SZL)  .  In  the  Zero  Equation  Model,  the 
definition  of  a  turbulent  viscosity  does  not  introduce 
another  differential  equation,  and  ANSYS  calculates  a 
characterist ic  length  applied  to  the  entire  model.  For 
this  model  ANSYS  applied  a  default  characteristic  length  of 
1  cm.  Neither  this  characteristic  length  nor  any  changes 
made  improved  upon  the  performance  of  the  default  SKE  model 
so  it  was  removed  from  further  consideration. 

The  GIR  model  was  more  stable  than  the  ZEM  model  but 
it  provided  results  significantly  different  than  the 
experimental  data.  It  relies  on  a  complex  formulation  for 
the  C„  term. 

The  SZL  model  is  a  simpler  model  than  NKE  or  GIR.  It 
estimates  a  low  level  of  turbulence  for  the  flowfield  and 
low  effective  velocity  that  can  quickly  become  unstable  in 
a  complex  flowfield.  The  C..  term  is  significantly  simpler. 

It  was  found  that  the  heat  transfer  did  not  steady  out 
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smoothly  and  it  consistently  overestimated  the  parameters 
used  to  determine  the  convection  coefficient,  h. 

Since  each  of  these  3  models  was  found  wanting 
compared  to  the  earlier  3,  they  were  not  pursued  in  much 
further  detail  in  this  study. 

F .  ADVECTION  SOLVERS 

For  this  simulation  the  ANSYS  recommendations  for 
solving  complex  flow  problems  were  used.  Single  solvers  by 
themselves  did  not  correctly  simulate  the  flows  but 
combinations  of  solvers  did.  For  the  first  30  iterations 
the  Monotone  Streamline  Upwind  (MSU)  advection  solver  was 
used  for  momentum,  turbulence,  pressure,  and  temperature. 
MSU  produces  first  order  accurate  solutions  and  is 
generally  very  robust.  In  the  current  study  in  particular, 
after  the  first  30  iterations,  density  was  permitted  to 
vary  and  the  MSU  was  still  used  for  momentum,  turbulence, 
pressure,  and  temperature.  For  iterations  60  through  90 
the  momentum  was  judged  stable  enough  to  switch  to  the 
Collocated  Galerkin  (COLG)  method,  an  exact  advection 
solver.  For  the  iterations  beyond  90,  the  COLG  solver  was 
used  in  the  solution  of  momentum,  turbulence,  pressure,  and 
temperature  quantities  for  best  consistency  and  stability 
in  the  final  results. 
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IV.  RESULTS  AND  DISCUSSION 


A.  NUMERICAL  EXPERIMENTS  CONDUCTED 

In  order  for  a  typical  numerical  run  to  be  considered 
successful,  which  could  then  be  compared  to  other  runs  with 
different  parameters,  a  number  of  numerical  performance 
criteria  had  to  be  met.  The  first  and  most  important  was 
grid  independence.  To  achieve  this,  the  number  of  elements 
in  the  model  was  increased  (in  steps  of  at  least  24000 

elements)  until  a  less  than  two  percent  change  was  observed 
in  the  total  heat  transfer  rate  between  two  consecutive 

runs.  Only  one  Reynolds  number  resulted  in  the  heat 

transferred  being  outside  the  2%  limit,  and  that  too  was 
less  than  3%. 

A  second  criterion  was  that  no  nodal  temperature  ever 
went  above  306. 4K  in  the  last  iteration  in  order  to  be 
consistent  with  the  imposed  boundary  conditions.  A  third 
was  that  the  ANSYS  generated  output  bulk  temp  was  within 

0 .  IK  of  the  calculated  temperature  based  on  the  mass  flow 
rate.  The  final  criterion  that  had  to  be  met  was  that  the 
sum  of  the  differences  in  nodal  temperatures  between  two 
consecutive  iterations  divided  by  the  sum  of  the  final 
nodal  temperatures  be  less  than  lxlCT4  to  ensure 
satisfactory  convergence. 

Only  when  all  these  criteria  were  met  was  a  run  deemed 
complete,  and  hence  ready  for  comparison  with  other 
completed  runs  . 

B.  TURBULENCE  MODEL  PERFORMANCE 

Three  models  (SKE,  RNG,  NKE)  that  were  able  to  meet 
the  above  robust  benchmark  criteria  were  compared.  As 
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noted  earlier,  three  others  (ZEM,  GIR,  SZL)  were  also 
considered  initially  but  failed  to  meet  the  basic  benchmark 
criteria  and  were  hence  not  pursued  further. 

1.  The  SKE  Model 

SKE  is  the  default  turbulence  model  implemented  in 
ANSYS.  The  model  generated  stable  solutions  and  in  order 
to  compare  it  with  Metzger's  experimental  work  an  array 
averaged  Nusselt  number  was  computed  for  each  variation  of 
each  model,  and  compared  with  its  corresponding  Metzger 
setup.  Baseline  results  using  SKE  as  the  turbulence  model 
consistently  overestimate  the  Metzger  Nusselt  number  with 
the  exception  of  the  high  Reynolds  number  run  in  the 
tightest  (X/D=1.5)  configuration  considered.  Various 
coefficients  in  the  SKE  turbulent  kinetic  energy  and 
dissipation  rate  equations  ( 1 . 1- . 2 ) could  be  varied.  Four 
out  of  a  possible  nine  coefficients  were  chosen  as  the  ones 
effecting  the  most  critical  changes  while  the  others  either 
dealt  with  buoyancy  terms  or  were  already  accounted  for  by 
varying  the  chosen  four.  Each  parameter  was  varied  by  an 
arbitrarily  chosen  value  of  30%  up  and/or  down  from  its 
default  value  and  a  run  was  made  with  only  that  parameter 
changed,  i.e.  all  else  was  held  fixed.  After  every  run 
benchmarks  were  checked  and  a  comparison  with  Metzger's 
Nusselt  number  was  made.  Table  4  summarizes  SKE 
performance  with 


respect  to  changes  in  Cl£  . 


Reynolds 

Number 

X/D 

S/D 

Variable 

Changed 

Percent 

Changed 

Simulated 

Nu 

Experimental 

Nu 

Percent 

Difference 

3980 

1.5 

2.5 

Default 

0 

33.84 

32.28 

4.83 

3980 

1.5 

2.5 

Cl-epsilon 

+30 

32.19 

32.28 

-0.28 

3980 

1.5 

2.5 

Cl-epsilon 

-30 

39.83 

32.28 

23.39 

7310 

1.5 

2.5 

Default 

0 

50.58 

49.61 

1.96 

7310 

1.5 

2.5 

Cl-epsilon 

+30 

50.67 

49.61 

2.14 

7310 

1.5 

2.5 

Cl-epsilon 

-30 

57.66 

49.61 

16.23 

13800 

1.5 

2.5 

Default 

0 

71.71 

77.75 

-7.77 

13800 

1.5 

2.5 

Cl-epsilon 

+30 

74.11 

77.75 

-4.68 

13800 

1.5 

2.5 

Cl-epsilon 

-30 

78.35 

77.75 

0.77 

3590 

2.5 

2.5 

Default 

0 

29.92 

26.73 

11.93 

3590 

2.5 

2.5 

Cl-epsilon 

+30 

24.58 

26.73 

-8.04 

3590 

2.5 

2.5 

Cl-epsilon 

-30 

36.13 

26.73 

35.17 

7340 

2.5 

2.5 

Default 

0 

45.69 

44.99 

1.56 

7340 

2.5 

2.5 

Cl-epsilon 

+30 

40.02 

44.99 

-11.05 

7340 

2.5 

2.5 

Cl-epsilon 

-30 

52.65 

44.99 

17.03 

13900 

2.5 

2.5 

Default 

0 

70.60 

71.61 

-1.41 

13900 

2.5 

2.5 

Cl-epsilon 

+30 

61.65 

71.61 

-13.91 

13900 

2.5 

2.5 

Cl-epsilon 

-30 

76.49 

71.61 

6.81 

Table  4 . 


Summary  of  SKE  performance  when  C{  is 
varied 


The  C]e  constant  precedes  the  viscous  dissipation  term 

in  the  dissipation  rate  equation  (1.2) .  Increasing  the 
coefficient  generally  results  in  a  lower  approximation  of 
the  Nusselt  number  than  the  default  case.  Specifically  the 
increase  in  C|E  contributes  to  an  increase  in  the  source 

term  of  the  dissipation  rate  transport  equation  and  results 
in  a  decrease  in  the  turbulent  kinetic  energy.  The  case  of 
Reynolds  number  of  13800  does  not  follow  the  same  trend  as 
the  other  Reynolds  numbers.  Figure  6  shows  this  decrease 
in  turbulent  kinetic  energy  as  displayed  by  ANSYS  6.0,  the 
view  is  from  the  adiabatic  plane  looking  down  towards  the 
lower  isothermal  plane.  The  top  figure  is  the  default  case 
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and  the  lower  figure  is  with  the  increased  ClE  multiplier 
(note  the  scale  of  the  colors) . 
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In  exploring  the  role  of  the  C2  coefficient,  it  must 
be  noted  that  C2  also  contributes  to  the  source  term  of  the 
dissipation  rate  transport  equation  (1.2)  in  a  nonlinear 
manner  and  it  was  found  that  a  decrease  in  C2  leads  to  a 
decrease  in  turbulent  kinetic  energy  at  each  node.  The 
trend  seen  in  Table  5  shows  that  decreasing  C1  leads  to 

lower  Nusselt  numbers  than  in  the  default  SKE  case,  again 
with  the  exception  of  the  case  of  Re  =  13800.  The  effect  of 
C,  on  the  model  increases  at  higher  Reynolds  numbers. 

Decreasing  the  turbulent  Schmidt  number  for  the 
transport  of  dissipation  rate,  Ot.  ,  affects  only  the 
diffusive  terms  of  the  dissipation  rate  transport  equation 
(1.2)  .  A  decrease  in  oe  leads  to  an  increase  in  the 
effective  diffusivity  of  the  dissipation  rate.  Table  5 
shows  that  a  change  in  Ot.  has  very  little  impact  on  the 
default  SKE  Nusselt  number.  The  model  at  Re=13900  was  run 
on  a  different  mesh  than  the  default  SKE  model  therefore  no 
conclusions  can  be  drawn  from  the  large  difference. 

The  turbulent  Schmidt  number  for  the  transport  of 
kinetic  energy,  o  ,  is  the  only  parameter  that  affects  the 

turbulent  kinetic  energy  transport  equation  (1.1)  through 
the  diffusive  terms.  A  decrease  in  o,  increases  the 
effective  diffusivity  of  the  turbulent  kinetic  energy. 
Table  5  shows  that  the  decrease  in  ak  has  a  very  small 

effect  on  the  default  Nusselt  number  similar  to  the  small 
effect  of  o  . 

e 
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Reynolds 

Number 

X/D 

S/D 

Variable 

Changed 

Percent 

Changed 

Simulated 

Nu 

Experimental 

Nu 

Percent 

Difference 

3980 

1.5 

2.5 

Default 

0 

33.84 

32.28 

4.83 

3980 

1.5 

2.5 

C2 

-30 

32.90 

32.28 

1.92 

3980 

1.5 

2.5 

Sigma-E 

-30 

33.86 

32.28 

4.89 

3980 

1.5 

2.5 

Sigma-K 

-30 

34.32 

32.28 

6.32 

7310 

1.5 

2.5 

Default 

0 

50.58 

49.61 

1.96 

7310 

1.5 

2.5 

C2 

-30 

49.83 

49.61 

0.44 

7310 

1.5 

2.5 

Sigma-E 

-30 

49.82 

49.61 

0.42 

7310 

1.5 

2.5 

Sigma-K 

-30 

49.95 

49.61 

0.69 

13800 

1.5 

2.5 

Default 

0 

71.71 

77.75 

-7.77 

13800 

1.5 

2.5 

C2 

-30 

73.59 

77.75 

-5.35 

13800 

1.5 

2.5 

Sigma-E 

-30 

71.29 

77.75 

-8.31 

13800 

1.5 

2.5 

Sigma-K 

-30 

71.17 

77.75 

-8.46 

3590 

2.5 

2.5 

Default 

0 

29.92 

26.73 

11.93 

3590 

2.5 

2.5 

C2 

-30 

23.89 

26.73 

-10.62 

3590 

2.5 

2.5 

Sigma-E 

-30 

29.06 

26.73 

8.72 

3590 

2.5 

2.5 

Sigma-K 

-30 

29.39 

26.73 

9.95 

7340 

2.5 

2.5 

Default 

0 

45.69 

44.99 

1.56 

7340 

2.5 

2.5 

C2 

-30 

40.00 

44.99 

-11.09 

7340 

2.5 

2.5 

Sigma-E 

-30 

44.88 

44.99 

-0.24 

7340 

2.5 

2.5 

Sigma-K 

-30 

44.88 

44.99 

-0.24 

13900 

2.5 

2.5 

Default 

0 

70.60 

71.61 

-1.41 

13900 

2.5 

2.5 

C2 

-30 

56.06 

71.61 

-21.71 

13900 

2.5 

2.5 

Sigma-E 

-30 

60.89 

71.61 

-14.97 

13900 

2.5 

2.5 

Sigma-K 

-30 

60.44 

71.61 

-15.60 

Table  5.  Summary  of  SKE  performance  when  C1 ,  Oe  , 

and  ck  are  decreased  by  30% 

For  the  SKE  model,  general  trends  observed  were  that 
for  the  lower  Reynolds  numbers  (3980,  7310,  3590,  7340)  the 

default  SKE  model  overestimated  the  Nusselt  number  for  the 
array.  For  higher  Reynolds  numbers  (13800,  13900)  the 
trend  was  reversed  and  the  default  SKE  model  consistently 
underestimated  the  experimental  Nusselt  number.  This  trend 
seems  to  be  independent  of  the  longitudinal  spacing  (X) 
between  pins . 

2.  The  RNG  Model 

The  next  model  analyzed  was  the  RNG  model.  For  each 
Reynolds  number  a  run  was  carried  out  for  the  default  model 

settings.  The  resulting  Nu  was  compared  to  the  same 
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experimental  Nu  from  the  SKE  runs.  There  are  10  model 
coefficients  in  the  RNG  model  that  can  be  modified  using 
ANSYS  6.0.  Only  two  of  these  parameters  (T)^  and  P  )  have  a 

different  effect  on  the  two  transport  equations  (1.1-2) 
compared  to  the  SKE  model.  Both  of  these  parameters  are 
the  multipliers  in  the  RNG  formulation  for  Cl  ,  equation 

(1.9) .  The  constant  T|m  has  the  greatest  effect  of  the  two 
on  the  Cj  formulation.  Hence  this  was  selected  as  the  only 
parameter  to  vary  in  the  RNG  trials;  the  others  were 
unchanged  from  the  SKE  model  since  they  had  been  accounted 
for  in  those  trial  runs.  The  value  of  in  the  RNG  model 

was  also  changed  by  30%  up  and  down  from  its  default 
setting.  Table  6  summarizes  the  resulting  changes. 


Reynolds 

Number 

X/D 

S/D 

Variable 

Changed 

Percent 

Changed 

Simulated 

Nu 

Experimental 

Nu 

Percent 

Difference 

3980 

1.5 

2.5 

Default 

0 

33.55 

32.28 

3.93 

3980 

1.5 

2.5 

Eta-inf 

+30 

36.80 

32.28 

14.00 

3980 

1.5 

2.5 

Eta-inf 

-30 

33.82 

32.28 

4.77 

7310 

1.5 

2.5 

Default 

0 

49.40 

49.61 

-0.42 

7310 

1.5 

2.5 

Eta-inf 

+30 

52.51 

49.61 

5.85 

7310 

1.5 

2.5 

Eta-inf 

-30 

50.88 

49.61 

2.56 

13800 

1.5 

2.5 

Default 

0 

71.00 

77.75 

-8.68 

13800 

1.5 

2.5 

Eta-inf 

+30 

73.33 

77.75 

-5.68 

13800 

1.5 

2.5 

Eta-inf 

-30 

75.13 

77.75 

-3.37 

3590 

2.5 

2.5 

Default 

0 

29.00 

26.73 

8.49 

3590 

2.5 

2.5 

Eta-inf 

+30 

32.63 

26.73 

22.07 

3590 

2.5 

2.5 

Eta-inf 

-30 

25.36 

26.73 

-5.13 

7340 

2.5 

2.5 

Default 

0 

44.25 

44.99 

-1.64 

7340 

2.5 

2.5 

Eta-inf 

+30 

48.45 

44.99 

7.69 

7340 

2.5 

2.5 

Eta-inf 

-30 

40.19 

44.99 

-10.67 

13900 

2.5 

2.5 

Default 

0 

68.32 

71.61 

-4.59 

13900 

2.5 

2.5 

Eta-inf 

+30 

72.58 

71.61 

1.35 

13900 

2.5 

2.5 

Eta-inf 

-30 

62.83 

71.61 

-12.26 

Table  6.  Summary  of  RNG  Performance  When  r)^  is 

Varied 
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The  RNG  results  for  the  default  model  underestimate 
the  experimental  Nusselt  for  all  but  the  lowest  2  Reynolds 
numbers.  This  trend  is  independent  of  the  X/D  ratio.  For 
the  lowest  Reynolds  numbers  the  RNG  default  model 
overestimates  the  Nusselt  number.  The  RNG  is  a  more 
sophisticated  model  in  that  it  takes  into  account  the 
variation  of  the  turbulent  kinetic  energy  in  the  flow  field 
when  determining  the  C1  term.  However  this  added  parameter 

provides  a  more  accurate  Nusselt  number  than  the  SKE  model 
only  at  Reynolds  numbers  less  than  or  equal  to  7340.  This 
agrees  with  prior  work  described  by  Shah  et  al  (2001)  for 
2-D  cases  using  the  RNG  model  in  similar  flow  scenarios. 
Varying  the  constant  does  not  have  predictable  results. 

Increasing  the  constant  by  30%  always  pushes  the  simulated 
Nusselt  number  higher.  Decreasing  the  constant  has  less 
predictable  results,  but  in  the  extended  configuration  (X/D 
=  2.5),  decreasing  T|to  decreased  the  Nusselt  number.  This 

appears  to  be  geometry  dependent.  The  results  of 
decreasing  can  be  seen  in  the  change  of  the  distribution 
of  turbulent  kinetic  energy  in  the  flow  field  for  the  3980 
Reynolds  number  in  Figure  7.  Figure  8  shows  the  difference 
in  the  temperature  field  for  the  same  conditions.  Both 
views  are  looking  down  through  the  adiabatic  plane. 
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Figure  7.  Change  in  the  turbulent  kinetic  energy 

distribution  of  an  RNG  run  at  Re=3980  where 
was  decreased  by  60% 
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Figure  8.  Change  in  the  temperature  distribution 

of  an  RNG  run  at  Re=3980  where  was  decreased 

by  60% 

In  other  runs  it  was  observed  that  any  change  up  or 
down  in  the  C2  constant  brought  a  change  in  the  Nusselt 

number  in  the  same  direction.  The  constant  C2  was  a 

consistent  predictor  but  its  effectiveness  decreased  the 
greater  the  change  in  the  constant.  It  was  unable  to 
change  the  Nusselt  numbers  by  enough  to  bring  them  closer 
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to  the  experimental  values  without  causing  the  transport 
equations  to  become  unstable. 

Other  runs  also  included  runs  where  C  was  decreased 

by  30%.  The  Nusselt  number  changed  in  the  same  direction 
for  all  Reynolds  numbers  for  a  given  spacing  configuration. 
When  X/D=1.5,  decreases  in  C„  brought  increases  in  the 

Nusselt  number  with  respect  to  the  default  RNG  runs. 
However  when  X/D=2.5,  decreases  in  C resulted  in  decreases 

in  the  Nusselt  number. 

3.  The  NKE  Model 

The  next  model  run  was  the  New  k  —  8  (NKE)  Model.  This 
model  is  proposed  in  the  literature  as  being  suitable  for 
high  Reynolds  numbers.  The  model  bears  this  out  for  its 
default  case.  The  difference  between  the  simulated  and 
experimental  Nusselt  numbers  decreases  as  Reynolds  number 
increases  but  only  for  the  cases  of  X/D=1.5.  At  Re=13900 
the  numerical  runs  were  unsuccessful.  The  results  are 
summarized  in  Table  7.  The  values  of  C,„  were  chosen  to  be 
evenly  spaced  between  0.43  and  1.0.  The  lower  number, 
0.43,  corresponds  with  the  default  setting,  and  1.0 

ri 

corresponds  with  the  limit  of  -  as  T)  goes  to  infinity. 

r|  +  5 
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Reynolds 

Number 

X/D 

S/D 

Variable 

Changed 

Value 

Simulated 

Nu 

Experimental 

Nu 

Percent 

Diff. 

MM 

IB 

2.5 

Default 

0.43 

35.94 

32.28 

11.34 

&SK 

m 

2.5 

Cl  M 

0.6 

31.44 

32.28 

-2.60 

3980 

1.5 

2.5 

Cl  M 

0.8 

32.88 

32.28 

1.86 

mm 

RB 

2.5 

Default 

0.43 

53.53 

49.61 

7.90 

H 1 

m 

2.5 

Cl  M 

0.6 

46.91 

49.61 

-5.44 

7310 

1.5 

2.5 

Cl  M 

0.8 

50.71 

49.61 

2.22 

13800 

1.5 

2.5 

Default 

0.43 

74.87 

77.75 

-3.70 

13800 

1.5 

2.5 

Cl  M 

0.6 

70.23 

77.75 

-9.67 

13800 

1.5 

2.5 

Cl  M 

0.8 

75.97 

77.75 

-2.29 

3590 

2.5 

2.5 

Default 

0.43 

31.66 

26.73 

18.44 

3590 

2.5 

2.5 

Cl  M 

0.6 

25.17 

26.73 

-5.84 

3590 

2.5 

2.5 

Cl  M 

0.8 

25.00 

26.73 

-6.47 

2.5 

2.5 

Default 

0.43 

46.45 

WESM 

3.25 

2.5 

2.5 

Cl  M 

0.6 

38.67 

-14.05 

7340 

2.5 

2.5 

Cl  M 

0.8 

39.07 

44.99 

-13.16 

13900 

2.5 

2.5 

Default 

0.43 

- 

13900 

2.5 

2.5 

Cl  M 

- 

13900 

2.5 

2.5 

Cl  M 

l—l 

- 

Table  7 .  Summary  of  NKE  Performance 

Changing  the  Cm  constant  for  the  X/D  case  of  1.5 

results  in  improved  performance  between  0.6  and  0.8, 
increasing  to  0.8  for  the  highest  Reynolds  number  of  the 
X/D  =  1.5  case.  In  the  X/D  =  2.5  case  changing  the  C1M 

constant  has  a  mixed  effect. 
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V. 


CONCLUSIONS  AND  RECOMENDAT I ON S 


A.  OPTIMIZATION  OF  THE  SHORT  PIN  FIN  HEAT  EXCHANGER 

Table  8  shows  the  optimal  constants  for  the  short  pin 
fin  heat  exchanger.  There  is  substantial  improvement  in 
the  X/D=1.5  spacing  ratio  model.  The  lower  Reynolds  number 
flows,  regardless  of  spacing  ratio,  are  not  well 
represented  by  the  SKE  default  model.  The  SKE  default  model 
consistently  overestimates  the  amount  of  turbulent  kinetic 
energy  in  low  Reynolds  number  flows  resulting  in  simulated 
flows  with  high  Nusselt  numbers.  The  increase  in  the  ClE 
constant  and  the  reduction  in  the  C„  constant  weight  terms 

in  their  specific  models  to  decrease  turbulent  kinetic 
energy  and  increase  the  dissipation  rate.  This  can  also  be 
seen  in  Table  4  where  Cle  is  increased  by  30%  and  this 
drives  the  dissipation  rate  up  as  the  transport  equations 
(1.1— .2)  are  solved  resulting  in  a  lower  Nusselt  number, 
closer  to  the  experimental  value. 


Reynolds 

Model 

Variable 

Combination 

Simulated 

Nu 

Experimental 

Nu 

SKE 

Def.  Nu 

Percent 

Improvement 

3980 

SKE 

Cl +30 

32.19 

32.28 

34.81 

96.44 

7310 

RNG 

Default 

49.40 

49.61 

50.58 

78.35 

13800 

SKE 

Cl -30 

78.35 

77.75 

71.71 

90.07 

3590 

RNG 

Cmu-30 

26.37 

26.73 

29.92 

88.71 

7340 

SKE 

Default 

45.69 

44.99 

45.69 

0.00 

13900 

RNG 

Eta+30 

72.58 

71.61 

70.60 

3.96 

Table  8.  Optimal  Constants  for  the  Cross  Flow 

Pin  Fin  Heat  Exchanger. 
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Table  9  shows  the  performance  of  the  optimal  constants 


against  Metzger's  values. 


Reynolds 

Number 

Model 

Variable 

Combination 

Simulated 

Nu 

Experimental 

Nu 

Percent 

Difference 

3980 

SKE 

Cl +30 

32.19 

32.28 

-0.28 

7310 

RNG 

Default 

49.40 

49.61 

-0.42 

13800 

SKE 

Cl -30 

78.35 

77.75 

0.77 

3590 

RNG 

Cmu-30 

26.37 

26.73 

-1.35 

7340 

SKE 

Default 

45.69 

44.99 

1.56 

13900 

RNG 

Eta+30 

72.58 

71.61 

1.35 

Table  9.  Optimal  Constants  Compared  with 

Metzger's  Experimental  Values 

B .  FUTURE  PLANS 

The  next  steps  for  this  study  should  be  to  explore 
changes  to  terms  in  the  models  that  affect  the  balance  of 
turbulent  kinetic  energy  and  the  turbulent  kinetic  energy 
dissipation  rate  (1.1-.2)  .  This  balance  seems  to  have  the 
most  significant  effect  on  the  performance  of  a  turbulence 
model.  Specific  changes  need  to  be  in  exploring  a  reduction 
of  the  Cj  constant  for  higher  Reynolds  numbers  in  the 

expanded  configuration  (X/D=2.5)  of  the  SKE  model.  There  is 
a  possibility  that  the  viscous  dissipation  term  becomes 
greater  as  the  flow  moves  farther  downstream  and  away  from 
turbulence  creating  obstructions  such  as  the  pins.  Table  4 
shows  that  the  changes  of  +/-  30%  of  C|e  bracket  the 

experimental  Nusselt  number  with  the  exception  of  the 

Re=7310  run.  This  leads  to  varying  constants  in  the  model 

that  are  directly  tied  to  source  and  dissipation  terms  as 

Cle  is.  Iterating  the  percent  change  should  bring  that 

constant  into  close  agreement  with  the  experimental  value. 

Another  specific  area  is  the  higher  Reynolds  number 
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performance  of  the  RNG  model.  The  constant  does  not 

provide  predictable  improvements  for  the  model.  At  higher 
Reynolds  numbers  there  may  be  a  benefit  in  increasing  the 
C2  term  that  reduces  the  contribution  of  the  viscous 

dissipation  rate,  which  appears  as  E2 . 

For  future  considerations,  temperature  differences 
must  be  increased  to  more  accurately  model  environmental 
temperatures,  H/D  ratios  should  be  varied  slightly  (0.5- 
2.0)  to  see  where  these  relations  need  to  be  improved. 

Another  environmental  concern  is  to  model  the 
rotational  environment  experienced  in  a  gas  turbine's  power 
turbine.  For  further  heat  transfer  performance  refinement, 
attention  will  have  to  be  paid  to  the  variations  available 
for  the  wall  turbulence  model. 

Many  compact  heat  exchangers  use  different  pin  shapes 
than  cylinders.  The  ANSYS  code  has  been  written  with 
vertical  planes  midway  between  each  row  that  will 
facilitate  changing  entire  rows  to  different  span  wise 
spacing  and  pin  shapes. 
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